Deep‐sea meiofaunal communities in the south‐eastern Levantine basin and their shaping factors – Morphological‐taxonomy‐free metabarcoding approach

Abstract The <3% dissimilar Amplicon Sequence Variant (ASV) clusters of the 18S‐V4 barcode were used as species‐proxies for the evaluation of ASV composition and ASV diversity indices characterizing the hitherto poorly investigated meiofaunal communities of the south‐eastern part of the Levantine basin. Accompanied by abundance measurements, the relationships of these characteristics with sedimentary and bottom terrain parameters were interpreted. The construction of community composition profiles, namely ASVs' list and their estimated abundances, was done using our previously established procedure (Harbuzov et al., 2022, Marine Genomics 65, 100980), combining metabarcoding with sample reads normalization by the abundance of hard‐bodied meiofaunal taxa. The study province included the 54–1418 m depth range, across vertical sub‐bottom horizons ranging 0–17 cm. Oxygen, hydrogen sulfide and methane concentrations in the pore water, as well as sediment grain size spectra and sedimentary protein and carbohydrate levels, were measured, followed by an evaluation of their involvement in the shaping of the meiofaunal communities' characteristics. Community composition was generally site‐and‐horizon dependent and its abundance decreased with increasing bottom depth and across sub‐bottom horizons, typical to benthic habitats which are nourished by organic carbon from the euphotic zone. The relatively sharply inclined continental slope bottom located in the northern part of the Israeli coast was an exception. Its meiofaunal community characteristics were speculated to be affected by intensive sediment mixing and lateral transport of food from the shelf, in addition to the effect of the euphotic zone‐originated food sources.


| INTRODUC TI ON
The present investigation was aimed at poorly studied aspects of soft substrate deep-sea habitat of the south-eastern Levantine basin but contains also methodological and ecological aspects of broader interest.Three issues were involved: (1) characterization of the meiofaunal communities abundances, species compositions and species diversities of the south-eastern Levantine basin, applying Amplicon Sequence Variants (ASVs) as species-proxies, (2) measurement of several sedimentary abiotic parameters across sub-bottom layers and evaluating their contribution to the shaping of the meiofaunal community characteristics and (3) further elaboration and implementation of novel molecular taxonomy analytical approach, composed of already known components, for characterizing the meiofaunal communities, compatible with species-level identification but independent of the difficult-to-perform morphological taxonomy (Harbuzov et al., 2022).The underlying hypothesis claimed that the biotic community characteristics at the studied sites are affected mainly by the availability of organic matter originating from the euphotic zone (Turner, 2015) with no other sources such as lateral food transportation or sub-bottom methane seepage (Smith, 2012).
It is widely accepted that biotic community studies require species-level identification, assuming that the species is the basic ecological biotic entity, uniformly affected by its surrounding conditions.Abebe et al. (2011) summarized the emerging difficulties while morphologically identifying nematodes, assumedly magnified if other meiofaunal groups were added.These difficulties included a global shortage of identification expertise and the practical need to boost the generally slow identification procedure while accomplishing large-scale ecological studies and environmental monitoring of meiofaunal communities.
Metabarcoding, an identification methodology based on DNA sequence taxonomy markers (termed barcodes) and High Throughput Sequencing (HTS) of DNA barcode reads in environmental samples, has been increasingly introduced to provide solutions to the meiofaunal identification difficulties.The various parameters that have to be considered while planning meiofaunal metabarcoding analysis were thoroughly discussed recently (Alberdi et al., 2018;Bruce et al., 2021;Gielings et al., 2021;Harbuzov et al., 2022;Pawlowski et al., 2018Pawlowski et al., , 2022;; Van der Loos & Nijland, 2020 and literature therein).Briefly, they include the choice of DNA source -whole sediment or sorted individuals; DNA extraction method; selection of barcode in relation to its variability in the studied taxa; adequate universal PCR primer pairs which could amplify the barcode from a variety of sampled taxa, partial or complete dependency of the reference library on morphological identification and last, sequence sources, regional or global.Harbuzov et al. (2022) presented and discussed the suitability of our metabarcoding approach and methodological choices in view of the practical and substantial challenges of the eastern Mediterranean and only recent modifications were detailed here.Our approach took advantage of the ability of the machine learning-based software DADA2 (Callahan et al., 2016) to distinguish unique species-compatible ASVs from Operational Taxonomy Unit (OTU) assemblies (Callahan et al., 2017).Therefore, our reference library was morphological taxonomy-independent.It was indicated by Harbuzov et al. (2022) that the number of ASVs in the reference library is roughly compatible with the number of inhabiting species.OTUs assemblies were HTSed only from DNA extracts of meiofauna collected in the studied province, assumed to create a comprehensive assembly containing all the local PCRamplified ASVs and almost all the inhabiting meiofaunal species.

The various soft bottom geomorphological features in the
Israeli part of the eastern Mediterranean were described in Kanari et al. (2020) and the literature therein and the bottom bathymetry is presented in Figure 1 including designation of major geomorphological features.The continental slope is characterized by canyons in the northern part of the coast (Almagor & Hall, 1984) and submarine landslides in its southern part (Katz et al., 2015).Two large slump complexes, Dor slump and Palmahim disturbance (Almagor & Garfunkel, 1979) were located on the slope.Five main morphologies were distinguished in the bathyal basin: folds, faults, sediment waves, deep-water channels and sediment fan lobes, formed by sediment transport processes and salt tectonics.Several regions of potential and actual methane seepage and brine drainage from sub-bottom layers were located in the continental slope and its borders with the bathyal basin (Eruteya et al., 2018;Herut et al., 2022;Lawal et al., 2022;Marine Ventures International, 2018).Last, the eastern Mediterranean is an ultra-oligotrophic sea (Berman-Frank & Rahav, 2012;Moore et al., 2013;Siokou-Frangou et al., 2010).

| Sampling, sorting, counting and abundance evaluation
Twelve sites situated within two transects roughly oriented perpendicular-to-the-shoreline were sampled in October 2018, one in the northern part of the Israeli coast, off Haifa and the second further south, off Tel Aviv, detailed in Figure 1.Canyon crests were sampled in the northern slope sites which were located in the canyon-rich region.
Additional Sites, HS1900, HS1822 and TA45 were used only for oxygen level evaluations, broadening our understanding of the sedimentary oxygen concentration patterns in the eastern Levantine basin.Three consecutive samples were collected in most sites, besides HS1418, in which only one sample was collected due to bad weather during the sampling cruise.Meiofauna was sub-sampled from a 0.25 m 2 box corer (BX-650, Ocean instruments, San Diego, CA) using a 9.4 cm diameter plexiglass core pushed down to the 17 cm horizon of the sediment.
Each core was horizontally sliced into 1, 2, 2, 2, 3, 3, 4 cm slices from the bottom surface downward, aimed at compensating for the decreasing number of individuals in deeper sub-bottom horizons.Processing of the sampled slices was carried out according to Harbuzov et al. (2022), ending with a mixture of isolated, intact and broken hard-bodied individuals and residual abiotic debris.Using Bogorov counting chambers under a stereoscope, the individuals in each sample were counted, sorted into the following major groups: Nematoda, Copepoda, Polychaeta, Isopoda, Ostracoda, Cumacea and Mollusca and abundances were calculated as a number of individuals/square meter.

F I G U R E 1
Sampling sites map.HS -Haifa transect; TA -, Tel Aviv transect, numbers on the right side of the site names designate depth in meters.TA45, HS1822, and HS900 were sampled only for oxygen evaluation.

| Grain size
Grain size spectra of the various slices in all sampled sites were carried out, one analysis per slice, by the Mastersizer 3000 (Malvern Panalytical) and were assembled into clay, silt and sand fractions, presented as Shepard's triangle (Shepard, 1954).The clay was determined as <4 μm particles, the silt between 4 and 63 μm and the sand as >63 μm particles.

| Protein concentration
Total sediment protein levels were evaluated according to Danovaro (2010), based on the colorimetric reaction of proteins with rameic tartrate and the Folin-Ciocalteau reagent in basic environment (pH 10).The procedure is detailed in Appendix 1.

| Carbohydrate concentration
Total sediment carbohydrate levels were evaluated according to Danovaro (2010), based on the colorimetric reaction between sugars and phenol in the presence of concentrated sulfuric acid (Dubois et al., 1956), optimized for sediments by Gerchakov and Hatcher (1972).The method is nonspecific and allows concentrations of total carbohydrates, cellulose included, to be determined.The procedure is detailed in Appendix 1.

| Sediment oxygen and hydrogen sulfide levels
A specially constructed device made of plexiglass by an engineering company (Aquazone, Israel) according to our specifications (Figure 2) was used to enable core-side measurements of oxygen and hydrogen sulfide levels and additional future parameters across sediment horizons down to 16.5 cm sub-bottom level.A sub-sample was cut from the box corer sample immediately on board using a 6 cm removable plexiglass core (Figure 2).Oxygen levels were evaluated using oxygen optodot sensors (PyroScience, GmbH, Germany) through the 3 mm wall of the removable core according to the manufacturer's instructions.Hydrogen sulfide levels were measured through holes in the core with a specific electrode (Unisense A/S, Denmark; type1 Sulf electrode) according to the manufacturer's instructions.

| Methane evaluations
Sediment samples (~2 mL) for methane concentration evaluations were collected with an edge-cut syringe from a plexiglass core perforated with 1 cm side holes at different horizons and immediately The device was used to measure abiotic parameters across a 16.5 cm vertical sediment core.The measurement instruments were connected to the power supply and the Manufacturers' computer programs.A core height adjustment wheel and a plastic piston were used to adjust core levels to the optodots and electrodes.
transferred into a glass bottle containing 5 mL of 1.5 N sodium hydroxide solution using gas chromatography according to Nusslein et al. (2003).

| DNA extraction
DNA from whole-sample sorted individuals was extracted using the E.Z.N.A.® Mollusk DNA Kit (Omega bio-tek, Cat-D3373), selected due to the presence of the cationic detergent cetyl trimethyl ammonium bromide (CTAB) in its lysis buffer, which improves DNA extraction from invertebrate tissues by efficient removal of mucopolysaccharides.Homogenization of samples in the lysis buffer was performed by the FastPrep bead homogenizer and its lysing matrix A beads (MP Biomedicals).It has to be noted that the pre-homogenized samples included also light-specific-gravity debris floated with the meiofaunal individuals during the density gradient centrifugation separation stage (Harbuzov et al., 2022).DNA levels were evaluated by fluorometry (QFX fluorometer, Denovix).
Unfortunately, the DNA preparations of the Tel Aviv transect slices were of low quality and many of them were not useful, although they served to establish a better DNA extraction method for the Haifa transect ones.The determination of the samples' ASV profiles included the below detailed steps.

| PCR amplification of sample's barcode assembly
The utilized barcode was region V4 of the small subunit rDNA, termed 18S.A PCR primer-pair, which produced a 470-490 bp PCR product was used by Harbuzov et al. (2022).High Throughput Sequencing (HTS) was performed from these PCR products applying the Illumina 300 × 2 bp platform enabling the sequencing of millions of OTUs belonging to meiofauna of the eastern Levantine basin (NCBI Sequence Read Archive (SRA) BioProject PRJNA791542).However, HTS done by several service laboratories using that Illumina platform elucidated relatively low high-quality yield due to the fast deterioration of quality score values toward the 3′ end of the 300 bp unpaired sequences, but the elucidated high-quality reads were used to design new primer-pair, producing shorter (325-445 bp) 18S-V4 barcode, tested for its suitability to be used to amplify barcodes from local meiofaunal samples following Harbuzov et al. (2022).The new primer-pair sequences were: Forward-5′-CCGCGGTAATWCCAGCHY-3′ and Reverse-5′-TTGGCAAATGCYTTCGCAKTHG-3′ coupled to the following Illumina overhangs: Forward-5′-ACACTCTTTCCCTACACGA CGCTCTTCCGATCT-3′ and Reverse-5′-GTGACTGGAGTTCAGACG TGTGCTCTTCCGATCT-3′.Each of the samples' DNA extracts were used as PCR template applying the Phanta Flash kit (Vazyme, China; cat.P520-00-AA) and a PCR design of 1′-98°C, (10″-98°C; 10″-58°C; 10″-72°C) X 35, 1′-72°C following the manufacturer's instructions.
PCR products that elucidated sufficient DNA levels with adequate fragment size were cleaned up using the NucleoSpin™ Gel and PCR Clean-up Kit (Macherey-Nagel, Germany, cat. No. 11992242) and were sent for HTS of ~150,000 reads per sample in a service laboratory (Syntezza Bioscience, Israel) using the Miseq device and the V2 2x250 cartridge (Illumina, USA).

| Metabarcoding analysis
According to our metabarcoding approach, a reference library has to be prepared from HTSs of local samples.Therefore, our existing reference library (Harbuzov et al., 2022) was combined with new local HTS results to construct an improved and more comprehensive 18S-V4 library, constructed according to the Harbuzov et al. (2022) procedure.Briefly, the FASTQ files of the present Haifa transect (NCBI bioproject PRJNA983550) and also HTS products resulting from smaller shelf projects, still unpublished, were filtered using the CUTADAPT software (Martin, 2011), including truncation of the poor 3′ sides of each sequence, primer removal and eliminating both short sequences (<180 bps) and sequences with >3 Ns in a row.Several quality scores between q = 20-35 were applied as CUTADAPT parameters and q = 25 was selected as the maximal one that still revealed a reasonable number of reads/slices.These filtered assemblies served as an input for the two components of the metabarcoding process.
The first component, the DADA2 analytical process (Callahan et al., 2016), applied through the Qiime2 software plugin (Bolyen et al., 2019)  were combined by repeating the clustering process, starting with the construction of a P-distance resemblance dissimilarity matrix using the Geneious Prime software (Biomatters LTD.) through its "tree" formation function, served in turn as input for clustering process of the sequences using the PRIMER-v7 software (Clarke et al., 2014;Clarke & Gorley, 2015) through its group average clustering protocol.Each ASV cluster with dissimilarity <3% among its member ASVs was considered by us a unique species-compatible entity and will be designated hereafter ASV-3.The reference barcodes were annotated by BLASTN against the nucleotide NCBI standalone database and non-18S sequences were eliminated (BioProject PRJNA983550).
The second metabarcoding component, pairing the CUTADAPTfiltered OTU reads of each of the Haifa transect slice samples was performed by the VSEARCH software (Rognes et al., 2016), using the VSEARCH merge-pairs plugin of Qiime2 (Bolyen et al., 2019).
Table 1 presents the average of the produced read numbers across the 83 successfully HTSed samples throughout the bioinformatic analytical process applying q = 25, showing the gradual reduction of read numbers throughout the metabarcoding process.
Metabarcoding was performed by the Geneious Prime software through its menus: 'Elign/assemble -map to reference', using 3% allowed read dissimilarity and 1% gaps of a maximum of 3 bp per gap, resulted in a profile of each slice, composed of the list of ASV-3s and the numbers of their mapped OTU paired reads.Only the ASV-3s annotated to meiofauna (Table 1) were further analyzed.The read numbers of each ASV-3 were normalized by the real number of counts in the sample, and a matrix of ASV-3s versus sampled slices was constructed.Each cell of the matrix contained an estimated count of a specific ASV-3 in each slice.The count normalization procedure is further detailed in the 'results' section below.

| Clustering of sample profiles
Count-normalized slice profiles were clustered by the PRIMER-v7 software (Clarke et al., 2014;Clarke & Gorley, 2015).The abundances were square root-transformed to reduce the effect of dominant taxa and the Bray-Curtis similarity measure was applied to the transformed data matrix.Slice clusters were visualized by applying the group average clustering method to create a dendrogram.The software PERMANOVA+ (Anderson et al., 2008) included in the PRIMER-v7 package was used for pairwise testing of the significance of differences among sample clusters under full model, applying type I sum of squares, a maximum of 10,000 permutations, unrestricted permutation of raw data and using the Monte Carlo correction.

| Calculation of ASV diversity indices
Diversity indices calculation was based on an entropy-like equation that included the Hill number variable (Chao et al., 2013(Chao et al., , 2014)).
Hill number = 0 evaluated only the number of ASV-3s while Hill number = 1 revealed the exponent of the Shannon-Weiner's Index which considered also ASV-3 abundances.The Chao et al. (2013Chao et al. ( , 2014) ) probabilistic approach enabled the calculation of observed and estimated ASV-3 diversity values, their standard errors, confidence limits, and their rarefaction and extrapolation curves, all used to estimate the diversity indices accuracy and the sampling sufficiency, by examining the proximity of observed and estimated diversity values from one index value per slice.Consequently, diversity indices were calculated from one combined sample profile per slice, constructed from three consecutive samples, making it more informative and accurate due to the larger sampled area probably covering more rare species.The combined abundances were normalized per 1 cm horizon and the calculation was accomplished by the R environment-based (R Core Team, 2021) software package iNEXT (Hsieh et al., 2016).

| Examining the correlation between biotic slice profiles and abiotic parameters
A matrix of the abiotic parameters sampled at each Haifa transect slice was constructed, including: total protein and carbohydrate concentrations, sediment silt/clay ratio, percentage of the sandy fraction of the sediment, oxygen concentrations in the pore waters water, depth and sub-bottom horizons.These variables were normalized prior to assessing their potential correlations with faunal composition (SR-transformed data) using the BIOENV procedure of PRIMER-v7 (Clarke & Ainsworth, 1993).The analysis made use of two resemblance matrices: the faunal composition matrix, calculated using the Bray-Curtis similarity index and the abiotic variables matrix, calculated using the Euclidean distance similarity index.The threshold of significant correlation was determined as the Spearman's correlation coefficient (rho) value which was higher than the maximal rho calculated from random permutations.

| General statistical methods
The R statistical environment (R Core Team, 2021) was used for general statistical tests.Statistically significant differences among groups of numerical values were tested by the pairwise Wilcoxon test or by the pairwise t-test and the Pearson equation was used for testing correlations.p < .05 was set as the threshold for null hypothesis rejection.

| Grain size
Grain size spectra of the various horizons in all sampled sites are presented in Figure 3 as Shepard's triangle, and for the most part, the results indicated a silt-clay bottom.In addition, the ratio of silt/ clay was calculated from each grain size spectrum (not shown).There were no significant differences of the silt/clay ratio among the Haifa transect samples and only two significant but minor differences in this ratio among the Tel Aviv ones.However, the average silt/clay ratio of the Haifa transect (1.08 ± 0.28) was significantly lower than that of the Tel Aviv one (1.67 ± 0.44), qualitatively observed also in the Shepard's triangle presentation, pointing out finer sediment off Haifa.The upper horizons of HS54 and HS771 revealed a substantial sandy fraction probably due to a high percentage of calcium carbonate shells (encircled in Figure 3).

TA B L E 1
Statistics of HTS read numbers throughout the analytical process.sites were not significantly different from each other.In the Tel Aviv transect, the protein levels of the shallower sites, TA76, TA155 and TA401 were significantly higher than those of the deeper TA941 and TA1173.TA1400 elucidated the significantly lowest protein concentrations in the transect.Generally, the deeper shelf and upper and mid slope on one hand and the deeper slope and bathyal basin sites on the other hand elucidated similar protein levels, the former were significantly higher than the later excluding the exceptional HS394.

| Carbohydrate concentration
Total carbohydrate concentrations in μg/g dry sediment were evaluated in the various sites and across the vertical horizons of both transects (Figure 5).Statistical tests of carbohydrate concentrations among sites and horizons revealed no difference across sub-bottom horizons in any of the sampled sites, and no discernible patterns were observed among both Haifa and Tel Aviv transects.Similarly, no significant difference was observed between the two transects.

| Oxygen levels
Oxygen levels in the pore water were evaluated in the various sites and across the vertical horizons of both transects (Figure 6).Three

| Hydrogen sulfide and methane levels
Methane seeps and high hydrogen sulfide levels were revealed at certain regions of the Israeli slope and its border with the bathyal F I G U R E 3 Shepard's triangle which represents the grain size characteristics along the Haifa and Tel Aviv transects and across horizons.The different horizons were not labeled for the sake of figure clarity.Outliers in the Haifa transect, which contain a higher sandy percentage are encircled.For site designation see Figure 1.

F I G U R E 4
Protein levels along the Haifa and Tel Aviv transects and across horizons.For site designation see Figure 1.
basin and more are predicted, as described in the Introduction section above.Therefore, the levels of both gases were evaluated in the studied transects.Residual average methane levels of 4.2 ± 4.8 μM were measured across all cores and no hydrogen sulfide residue was found in any of the sampled cores.

| Meiofaunal abundance
Table 2 presents the total numbers of higher-than-species counts of hard-bodied species, summed to 56,231 specimens, aimed at demonstrating the real number of handled individuals.Figure 7 presents the meiofaunal abundances of these individuals in the two transects and across the various horizons.Occasional broken individuals and the disintegrated soft-bodied taxa compromised the counting accuracy (see details below).Nevertheless, the counting provided a rough estimate of the meiofaunal abundance in the various slices.HS394 abundance was significantly higher than all other Haifa transect sites excluding HS54 which, in turn, was similar to all other sites excluding minor differences from HS1418.In the Tel Aviv section, TA76 abundance was significantly higher than all other sites besides TA155, which has in turn similar abundance with all F I G U R E 5 Carbohydrate levels along the Haifa and Tel Aviv transects and across horizons.For site designation see Figure 1.

F I G U R E 6
Oxygen level in sediment pore water along the Haifa and Tel Aviv transects and across horizons.For site designation see Figure 1.

TA B L E 2
Higher-than-species counted hard-bodied taxa summed across all replicate cores and horizons in each site.For site designation see Figure 1. the deeper sites.There was no significant overall average abundance difference between the two transects.Generally, excluding HS394, deeper shelf and upper slope abundances were higher than those of the deeper slope and bathyal basin.No individuals were elucidated at the deeper horizons of the Tel Aviv transect (Figure 7).

| Reference library
The reference library included 4154 ASVs, clustered into 2133 ASV-3 clusters.All the sequences were annotated and their families or above-family similarities were used for their designation.It has to be emphasized that the ASV designations were not considered exact identifications but a unique label for each ASV in our NCBI deposited library (BioProject PRJNA983550), which only indicated ASV taxonomy relatedness (e.g.Spionidae_32 resembles the sequence of a spionid polychaete in the NCBI collection and the number at the right side mark the unique ASV).The accuracy of an ASV designation is assumed to depend on the representation intensity of its taxon in the public databases.Hence, there could be less exact taxonomic assignments.The majority of the ASV annotations belonged to marine infaunal species, Protozoa and Metazoa.A number of 1353 ASV-3s were annotated to meiofaunal taxa, and were implemented as the reference library used for the metabarcoding analysis (Table 3).

| Sample and site ASV compositions
The resulting slice profiles were combined into an ASV × slice matrix which included the estimated counts of each ASV-3 in each slice.It Meiofaunal abundances in the various sampling sites and across horizons, in 1 m 2 area normalized to 1 cm horizon thickness.For site designation see Figure 1.The Y-axis scale is logarithmic (log 2 ).
TA B L E 3 Division of the ASV-3s according to their annotations and their number.

Functional groups Taxa
No has to be noted that the number of counted individuals was compromised due to broken individuals and to the relative fragility of the soft-bodied Platyhelminthes, Gastrotricha, Gnathostomulida and Xenacoelomorpha which were not counted at all, but comprised an average of 8.9% of the mapped reads.In spite of these potential counting inaccuracies, the normalization was postulated to better represent real ASV abundances than the relative abundanceindependent paired-OTU read numbers.
The (counts/number of ASV-3s) index, calculated for each slice profile was expected to be >1 if those reads would have originated only from counted individuals.Figure 8a demonstrated that the index is <1 in part of the normalized profiles, indicating background sources of ASVs not emerging from counted individuals but from compromised counting as described above and also from pray DNA in the digestive system of predators and DNA adsorbed to solid items of the sorted samples: low specific-gravity debris and fauna.Seawater-dissolved DNA is less likely to remain in the DNA extracts due to the preceding washing steps, but its presence could not be entirely ruled out.To partially correct this bias, an empirical semi-arbitrary approach was applied.The ASV-3s × slices matrix cells with low counts were gradually zeroed to a point where the (Counts/number of ASV-3) index was >1 for all the slices, revealing a range of 1-33 Counts/number of ASV-3 index across the entire matrix and a total of 990 ASV-3s (speciesproxies) across all sites.
Figure 9 presents the dendrogram constructed by cluster analysis of all the normalized and modified meiofaunal slice profiles.
The black dots in Figure 9 designate clusters that are significantly different from each other, examined by PERMANOVA+ (p < .05).
Generally, the profiles were assembled according to sampling sites and horizons, excluding the profiles of HS394 that were assembled together across horizons (Table 4).The outliers, designated by rectangles in Figure 9, were divided into two types: slice profiles that were feebly clustered to others (<20% similarity) and those that were out of their spatial context emphasizing slices from site HS771.
Part of the eight outliers in the dendrogram of Figure 9 may result from a relatively small number of ASVs, probably caused by insufficient sampled area, mostly from deeper sites or deeper horizons.No explanation could be provided for other outliers.Generally, diversity is decreasing in deeper sub-surface horizons.

| ASV diversity
However, the noticeable exception is the mid-slope site HS394 peaking at the 6 cm horizon.Lower peaks were observed for sites HS122, HS771 and HS1134 at the 2, 8.5 and 8.5 cm horizons, respectively.
All these peaks are weakly demonstrated also in the Exp Shannon- F I G U R E 9 Clustering of ASV meiofaunal profiles from the various slices.Outliers are framed by rectangles demonstrating clusters at a low similarity level (red rectangles) or slices that were clustered outside their spatial context (blue rectangles).Dashed vertical lines mark the 20% and 40% similarities.Black dots mark significantly different unique slice profile assemblies.The name of each slice is composed of the transect name (HS -Haifa transect), the sampling depth in meters (e.g., HS771), the consecutive sample at each site (e.g., HS771C), and short designation of horizon vertical depth in cm when 0=0-1, 1=1-3, 3=3-5, 5=5-7, 7=7-10, 10=10-13 and 13=13-17 cm.
HS394, HS771 and HS1134.Pairwise percentages of ASV-3s' sharing among sampling sites are presented in Table 5, showing that although meiofaunal communities were generally site-specific, around 18%-32% of the ASV-3s were present in more than one site.ASV-3s from site HS1418 were not included in the calculations of shared ASV-3s due to poor quantitative sampling.

| Correlation with abiotic conditions
Correlations of sample profiles with several abiotic parameters of the Haifa transect were tested by the BIOENV tool of PRIMER-v7, presented in Table 6.Only spatial parameters: sub-bottom horizons and bottom depths were significantly correlated with biotic slice profiles.

| Technical aspects
Three parameters were evaluated for the characterization of the deep-sea meiofaunal communities of the Israeli Mediterranean waters: abundance, ASV-3 composition and ASV-3 diversity indices.
Soft-bodied meiofaunal taxa were identified in the metabarcoding results but were absent in the sorted and counted individuals, probably due to the applied sorting methodology.A more adequate sorting methodology for these taxa was suggested by Balsamo et al. (2020) and it would have to be applied in the future for elucidation of that missing meiofaunal fraction.
Non-meiofaunal ASV-3s which were elucidated through their annotation were not the aim of this study and most of their related TA B L E 4 Division of slice profiles into clusters (see Figure 9).For site designation see Figure 1.
individuals were probably lost during the sorting process.However, many of their barcodes were elucidated in the reference library as a qualitative by-product (Table 2 and BioProject PRJNA983550).Their variety initially indicated the suitability of the utilized PCR primerpair for evaluation of their species composition upon the application of the more suitable DNA extraction method from unsorted sediment samples.

| Metabarcoding aspects
Metabarcoding introduces two theoretical inherent and incurable biases to the quantitative evaluation of species relative abundances: A disadvantage of complete detachment of the species identification from morphology is the loss of structure-related functional information.Assignment of morphological structures to ASVs could, at least partially, be restored by annotating the ASVs to public DNA sequence databases and when achievable, also the performance of morphological identification, not necessarily at the species level, similar to the abovementioned annotations, accompanied by amplification of the ASV from the identified individual.These higher-than-species identifications could provide reasonable structure information in most cases.
The selected 18S-V4 barcode elucidated a high level of ASV variability in the communities of the studied region, hence, in addition to the indicated species compatibility (Harbuzov et al., 2022) the used barcode provided a reasonably diverse list of species-like entities enabling the study of the meiofaunal community ecology.
The semi-arbitrary 3% dissimilarity threshold that was applied here to distinguish between species-compatible ASV clusters, fairly agrees with Wu et al. (2015) who studied 18S variability in copepods.Attempts to evaluate ranges of barcode variabilities including the one of 18S, mainly based on sequences downloaded from public databases, provided ranges and also revealed outlier species with exceptionally high variability (Phillips et al., 2019;Tanabe et al., 2016).
However, considering the complexity of testing barcode variability across a broad range of species, it is suggested to apply an empirical approach for the determination of appropriate cluster threshold across multi-species samples, by performing a metabarcoding procedure of template sample HTSs against alternatives of a reference library clustered using several threshold percentages and examining the number of resulted ASVs to recognize the patterns of changes in ASV numbers in relation to the threshold percentage.This analysis deserves a separate bio-informatic study.

| Ecological aspects
The continental slope of the eastern Mediterranean was recently observed to contain methane cold seeps and high levels of sulfides <0.21 at several, not entirely determined, areas (Eruteya et al., 2018;Herut et al., 2022;Lawal et al., 2022;Marine Ventures International, 2018).
The present study was aimed at soft bottom regions which were assumed to have a low probability of containing these emissions.
However, as a precautionary measure, methane and hydrogen sulfide levels were evaluated at all sampled horizons and indeed were found to be negligible.
Although the meiofaunal communities were significantly different among sites, they share a considerable amount of ASV-3s (Table 5 and text in the results section) and the significantly different site-related sample profiles were partially a result of different relative abundances rather than different ASV-3s.
The two diversity indices, ASV-3 richness and the exponent of the Shannon-Weiner of HS54, HS122 and HS1134 elucidated gradually downward decreased values till the 6 cm horizon from which their values remain constant (Figure 10).In view of the significant correlation of the two indices with abundance (Figure 7a) across horizons in these three sites it is speculated that relatively low abundance and consequent insufficient sampling caused the disappearance of rare ASVs.Decreased oxygen levels in deeper horizons (Figure 6) are less likely the reason for the decreased number of ASV-3s, as only the HS1134 levels were significantly correlated to the diversity indices.
The lack of correlation between the overall site ASV-3 richness and the exponent of the Shannon-Weiner index (Figure 11) may be explained by less even ASV-3 distribution patterns causing decreased exponent of the Shannon-Weiner index mainly in HS394.The complex behavior of the ASV diversity indices of HS394 and HS771 will be discussed below in a broader context.
Our working hypothesis which claimed the euphotic zoneoriginated food source as the main meiofaunal abundance shaping factor was generally true, indicated by the decrease of abundance with depth and across vertical horizons.ASV-3 compositions in the Haifa transect were related to depth as revealed by the BIOENV analysis, interpreted by us as an indirect depth effect, probably resulting from the decrease in available food sources with depth.However, this abundance decrease was not linearly related to depth but also to the bottom terrain and certain sites showed statistically similar abundances although sampled at different depths.This hypothesized terrain effect in the Haifa transect was based on the bottom steepness in the continental slope (Figure 1), the fluidic, semi-liquid mud visually observed in slope samples, emphasizing HS394, but also in HS771, the downslope currents that were speculated to flow in this canyon-rich region (Kanari et al., 2020) and the intermittent downslope turbidity currents that were observed there (Jaijel et al., 2023).These combined factors may cause sediment vertical mixing and downward sediment transport.Indeed, the more even meiofaunal distribution across horizons (Figure 7), the peaking ASV-3 diversity indices in deeper horizons (Figure 10) mainly in HS394 but also in HS771 and the HS394 slice profile cluster (Figure 9) which is indifferent to vertical horizons, indicate vertical mixing.The significantly high protein levels in HS394 (Figure 4) and also the generally higher average protein level in the Haifa transect than that of the Tel Aviv one may indicate organic matter transportation from the shelf and the upper slope, driven by downward currents and accompanied by significantly abundant meiofauna in the Haifa transect slope.
The carbohydrate and protein concentrations did not elucidate the same pattern in relation to depth and transect.Unlike the protein concentration, no trend could be determined in the carbohydrate levels and only a range of carbohydrate concentrations could be defined.It is hypothesized that proteins originated mainly from live organisms, whereas carbohydrates may originate from both live individuals and detritus.Total carbohydrate and protein levels were reported by Danovaro et al. (1999) from the Aegean Sea.The total carbohydrates were much higher in Danovaro et al. (1999) than the presented ones here.The existence of spatial differences in protein levels and the lack of carbohydrates spatial differences are similar in both studies, ours and Danovaro et al. (1999).
Pore-water oxygen levels in the studied sites elucidated similar characteristics to other soft substrates characterized by welloxygenated water above the bottom.Namely, fast biologically related depletion of oxygen toward deeper horizons in shallower waters is contrary to higher oxygen levels in deeper horizons of deeper sites (Glud, 2008).Although abundance decreased in deeper horizons, meiofaunal specimens were observed in apparently oxygendepleted ones, indicating adaptation of the meiofauna to low oxygen levels.

| Comparative ecological aspects
Several deep-sea studies of meiofaunal communities that were carried out in the eastern Mediterranean at relevant depths were compared to the present study (Table 7).The minimal size of sorted individuals slightly differed among the compared studies (20-45 μm) and similarly, the sampled range of vertical horizons (6-16 cm) which may slightly compromise the comparisons.Table 7 comparisons demonstrated lower abundance in the south-eastern edge of the Levantine basin than that of the western part of it, interpreted by the west-east primary production decrease across the Mediterranean (Siokou-Frangou et al., 2010).The numbers of nematode ASV-3s revealed during the present study showed the same order of magnitude in comparison to the morphological estimates done in other studies, further indicating the ASV-3s as species-proxies.

| Future recommended studies
During the sorting and the bioinformatic analyses presented here, several decisions were made, part of them were semi-arbitrary.The major processing improvement should be the elucidation of softbodied taxa.Several reasonable but semi-arbitrary decisions were made throughout the analytical procedure: applied quality score level, the threshold of counts/number of ASV-3s index and the threshold percentage of species-compatible ASV cluster.Their optimization efforts should be further advanced.
The meiofaunal communities of three more provinces of the Israeli Mediterranean deep-water soft substrate are required to be characterized for obtaining more accurate habitat delimitations required for environmental regulation of this intensively anthropogenically used region (Tom et al., 2023) using the already applied methodologies.These are the southern part of the slope (Tel Aviv transect), more moderately inclined than the northern transect, lacking canyons and containing lower sedimentary protein levels, the methane and sulfide seepage regions and the western deeper bathyal basin down to 2000 m containing considerable oxygen levels at deeper horizons (Figure 1).
An interesting and multi-disciplinary issue to be studied is the contribution of lateral transport and vertical sediment mixing processes on the slope meiofaunal community characteristics, involving the speculated effect of steep terrain and intermittent and permanent downslope currents.

| CON CLUS IONS
In summary, the study provided a novel detailed description of the meiofaunal community of a deep-water soft substrate in the Note: The table is arranged according to increasing depths and specific habitats (chemosynthetic sites and the Eratosthenes mountain) were avoided.
created a list of unique paired ASVs.Short sequences (<325 bp) were manually removed from this new reference library, realized to be erroneously paired or containing only partial sequences of the barcode.The already available ASVs of the Harbuzov et al. (2022) library were trimmed to fit the length and the genomic location of the new barcode and the two ASV libraries, old and new,

3. 2
| Protein concentrationTotal protein concentrations in mg/g dry sediment were evaluated in all the slices (Figure4).Statistical tests of protein concentrations among sites and horizons revealed no difference among sub-bottom horizons in any of the sampled sites.Protein concentrations were significantly higher in the Haifa transect than in the Tel Aviv one by an average of 1.43 mg/g dry sediment.Pairwise comparisons of sites in each of the transects revealed the following statistics: in the Haifa transect, HS394 protein concentrations were exceptionally high and different from all other sites.HS1418 concentrations were the lowest in the transect and significantly lower than all sites excluding HS1134.The protein levels of the HS54, HS122 and HS771 measurements were done for each site and horizon in the Tel Aviv transect, whereas only one or two were measured in Haifa transect due to technical failures.The oxygen levels were evaluated also in three sites outside our faunistically sampled ones and they are presented here as they expand the depth range of evaluated oxygen level patterns.A faster depletion of oxygen toward deeper horizons was observed in the Haifa transect in comparison to the Tel Aviv one and substantial oxygen levels were observed in the bathyal basin in both transects even at the 15.5 cm horizon of the deeper sites.
Estimated ASV-3 richness and the exponent of the Shannon-Weiner index were calculated for each Haifa transect site across the various horizons (Figure 10a,b).Standard error values of estimated indices were elucidated and rarefaction and extrapolation graphs with gray areas of 95% confidence limits were produced (Appendix 2), demonstrating sufficient sampling effort with almost identical observed and estimated ASV-3 diversity values, negligible standard errors, and a visible asymptote in all the rarefaction graphs with very narrow area of the 95% confidence limits.
Weiner index graphs.The two ASV-3 diversity index values across horizons were significantly correlated between themselves and also to their slice abundances, besides the indices of HS394 that were not correlated neither between themselves nor to the slice abundances.Only the HS1134 indices were significantly correlated to horizon-compatible oxygen levels.ASV-3 richness and the exp Shannon-Weiner indices of combined horizons at each site are presented in Figure 11, revealing a continental slope peak of species richness in the mid-slope HS394 site but a deep of the exp Shannon-Weiner index in the same site indicating uneven ASV-3 abundances in the HS394 profile with a community characterized by several dominant species.No significant correlation was shown between the two indices.A number of 121 ASV-3s were shared by deep shelf and slope sites and 76 ASV-3s were shared by the samples from HS54, HS122, F I G U R E 8 Counts versus Counts/ number of ASVs scatter graphs.Each black dot in the two graphs represents a slice.(a) Slices with high background reads, (b) Slices with reads resulted, at least partially, from whole individuals.
Misplaced slice profiles are designated in their cluster of presence.F I G U R E 1 0 Estimated ASV diversity indices in the various sites and across horizons.(a) ASV richness, (b) exponent of the Shannon-Weiner index.Estimated standard errors calculated by iNEXT software were very small and could not be observed on the graph, although introduced.HS1418 was poorly sampled and its results are probably an underestimate.For site designation see Figure1.ASV diversity indices per Haifa transect site.(a) ASV richness, (b) exponent of the Shannon-Weiner index.Estimated standard errors calculated by the iNEXT software were very small and could be barely observed only on graph b, although introduced.HS1418 was poorly sampled and its results are probably an underestimate.For site designation see Figure1.
potentially different PCR amplification efficiency for different species-specific DNAs due to different mismatch levels between the primer-pair and the DNA of various species and different copy numbers of the barcode in different species and individuals.Another and at least partially curable bias results from the assignment of relative read numbers to each ASV and the consequent lack of correlation with actual counts of individuals.The normalization of the read profiles by absolute counts was advisable to partially restore estimated counts and improve the accuracy of clustering among sample profiles.The subtle division among morphologically none identified ASVs, the soft-bodied individuals that were present in the metabarcoding profile but not in the physical counting, and also broken hard-bodied individuals, render this normalization less accurate, but despite the above-described inaccuracy sources, the metabarcoding-related evaluation of species composition presented here is a preferred practical way to accomplish community meiofaunal analyses, especially in faunistically unexplored habitats.On top of the practical necessity, HTS reads and ASVs are indifferent to taxonomy abilities, to disputes among experts, and to objective lack of taxonomy knowledge, rendering the comparability of sample profiles over time and space more accurate and repeatable in comparison to methods that rely partially or completely on morphological identification.This trait is especially important for environmental biotic monitoring or other time series studies.

HTSed pairs CUTADAPT -filtered pairs - Quality score = 25 Paired reads (VSEARCH) Reads mapped to the reference library Reads mapped to the meiofaunal part of the reference library
TA B L E 7 Comparisons of meiofaunal abundances and a number of species between the present article and others were performed in relevant depths of the eastern Mediterranean.